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COMPUTATIONAL  HYDRODYNAMICS: 
AN  INTERACTIVE  COMPUTER  CODE  FOR 
WATER  WAVES  IN  OPEN  CHANNELS 


I.  INTRODUCTION 

The  analysis  and  solution  of  free  surface  wave  problems  requires  the  development  of  a  model 
which  can  represent  the  physical  phenomenon  under  study.  However  the  complexity  of  the  equations 
describing  surface  waves  imposes  certain  restrictions  in  the  type  of  solution  that  one  can  use.  Most 
often  these  equations  are  solved  by  some  approximate  method  which  requires  the  development  of  a 
computational  model  implementing  such  a  method. 

For  the  purpose  of  studying  water  waves  in  open  channels  a  computational  model  has  been 
developed.  The  analytical  formulation  and  the  description  of  the  model  can  be  found  in  [1].  The  one¬ 
dimensional  form  of  the  shallow  water  equations  is  solved  numerically  by  the  finite  element  method 
(FEM).  The  simplest  form  of  approximation  in  space  and  time  is  utilized  for  the  FEM  and  the  derived 
model  is  capable  of  producing  accurate  solutions  to  a  variety  of  problems  with  free  surface  motion 
[2,31. 

The  numerical  solution  of  the  system  of  equations  which  represent  the  computational  model  is 
obtained  by  a  computer  code  developed  on  an  HP  1000  computer.  A  modification  of  the  code  was 
neccssaiy  in  order  to  utilize  the  code  in  an  interactive  way  and  to  be  able  to  solve  different  types  of 
cases  with  the  least  amount  of  user  effort.  Following  a  brief  description  of  the  analytical  model  in  the 
first  part  of  this  report,  the  computer  model  is  presented  in  the  second  part.  The  different  parts  of  the 
computer  model  are  given  in  a  flow  chart  diagram  and  for  each  part  of  the  required  subroutines  are 
described.  In  the  last  pan  of  the  report  the  solutions  of  some  example  problems  are  given.  A  user's 
manual  and  a  listing  of  the  computer  code  are  included  in  Appendices  A  and  B  respectively. 

II.  COMPUTATIONAL  MODEL 
1.  Analytical  Medal 


The  motion  of  a  body  of  water  with  a  free  surface  in  a  channel  can  be  described  by  a  set  of  two 
nonlinear  partial  differential  equations  (shallow  water  equations).  The  one  dimensional  form  of  these 
equations  is 
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where  h  is  the  instantaneous  water  surface  elevation  from  the  mean  surface  level,  H  is  the  total  depth 
of  the  fluid,  V  is  the  depth-averaged  velocity  in  the  x-direction.  r  is  the  time  and  g  is  the  gravitational 
acceleration.  Equations  (1)  and  (2)  are  cast  into  a  variational  form  and  the  FEM  is  applied  to  derive  a 
system  of  ordinary  differential  equations  given  by 
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where  (A/]  is  the  mass  matrix,  t£]  is  the  convective  matrix,  {Q}  is  a  vector  representing  boundary 
forces  and  {«)  is  the  vector  representing  either  {/»}  or  { V). 

The  system  of  equations  represented  by  Eq.  (3)  is  reduced  to  an  algebraic  system  by  approximat¬ 
ing  the  time  derivative  and  the  system  then  is  given  by 

+  (4) 

where  the  matrices  [A  1  and  (5]  are  combinations  of  the  matrices  [A/]  and  [AT].  The  vectors  {q"+l}  and 
{?”}  represent  the  unknowns  at  time  levels  (n  +  1)  and  (n)  respectively.  The  derivation  of  the  above 
equations  can  be  found  in  [1]. 

The  system  of  equations  represented  by  Eq.  (4)  is  solved  numerically  by  a  computer  code 
developed  at  NRL  for  that  purpose.  The  different  parts  of  the  code  are  described  in  the  following  sec¬ 
tion. 

2.  Cenpvter  Model 

The  analytical  model  described  in  the  previous  section  can  be  used  to  solve  a  variety  of  physical 
problems  concerned  with  water  waves.  Solutions  are  obtained  by  casting  the  analytical  model  into  a 
computer  model.  A  general  description  of  the  computer  model  is  given  in  the  flow  chart.  Fig.  I.  The 
different  stages  of  the  model  generate  the  required  information  and  perform  the  numerical  solution  for 
the  specified  problem. 

In  the  first  stage  the  physical  problem  under  consideration  is  defined.  This  includes  the  particular 
channel  geometry,  boundary  conditions  and  initialization  of  all  program  parameters.  The  subroutines 
called  in  this  stage  are  the  INPUT  and  GEOMT  subroutines.  In  the  next  stage  the  system  of  algebraic 
equations  is  formed  (subroutine  ELMAT)  and  the  system  is  solved  (subroutine  QSOLV)  for  the  partic¬ 
ular  type  of  boundary  conditions  under  consideration  (subroutines  MPART,  DCOMP,  and  SOLVE). 
Following  the  solution,  the  values  for  the  surface  elevation,  velocity  and  pressure  are  computed.  The 
process  is  repeated  for  every  time  increment.  The  last  stage  of  the  code  performs  such  tasks  as  output 
of  results,  plotting  of  results  either  at  specified  time  intervals  or  at  the  end  of  the  solution.  This  stage 
includes  subroutines  XOUT,  XPLOT,  and  TOUT,  TPLOT. 

As  was  mentioned  earlier  the  noninteractive  version  of  the  computer  code  can  be  used  on  any 
computer.  The  interactive  version  of  the  code  utilizes  the  HP- 1000  graphics  terminal  2648  and  certain 
changes  will  be  necessary  in  order  for  the  code  to  be  used  on  another  machine. 

The  code  has  been  optimized  for  minimum  storage  requirements  and  execution  time.  Further 
optimization  can  be  achieved  by  utilizing  the  vector  instruction  set  (VIS)  of  the  HP- 1000  computer  sys¬ 
tem.  Although  this  reduces  the  execution  time,  it  also  restricts  the  transferability  of  the  code  to  other 
machines. 

The  description  of  each  of  the  subroutines  is  given  in  the  next  section. 

III.  COMPUTER  PROGRAM 

3.  Pragnun  Description 

The  computer  program  is  composed  by  the  main  program  (program  KYMA)  and  eleven  subrou¬ 
tines.  The  dimensions  of  the  arrays  are  set  for  a  maximum  of  222  equations.  This  allows  for  III  nodal 
points  and  two  unknowns  per  point.  The  storage  requirements  for  the  set  dimensions  is  about  26  K 
words  or  26  pages  in  the  HP- 1000  system.  If  a  larger  number  of  nodal  points  is  required,  then  the 
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extended  memory  (EMA)  of  the  system  has  to  be  used  for  additional  storage.  The  operating  system  of 
the  HP- 1000  computer  on  which  the  code  is  running  is  the  RTE-IVB. 

PnpmKYMA 

This  is  the  main  part  of  the  code  and  is  used  as  a  driver  for  the  calling  sequence  to  the  different 
subroutines.  The  definitions  of  the  names  of  the  variables  and  parameters  used  in  the  code  are  given  in 
the  listing  of  the  code  in  the  appendix.  The  first  parameter  that  needs  to  be  defined  is  the  logical  unit 
(LU)  number  for  the  input  device,  i.e.,  LU  ™  l,  the  input  device  is  the  users  terminal.  The  call  to 
INPUT  sets  the  parameters  and  defines  the  problem  to  be  solved.  Next  the  user  has  to  specify  the  time 
for  starting  the  output  of  the  results  (TSTART)  and  the  time  to  stop  this  task  (TSTOP).  If  TSTART 
<  TIME  <  TSTOP  then  subroutines  XOUT  is  called.  For  each  time  increment  a  call  to  QSOLV  js 
required  for  the  solution  of  the  system  of  equations.  This  process  is  continued  until  the  solution 
reaches  the  specified  maximum  time  for  the  simulation. 

INPUT 

This  is  the  primary  parameter  definition  subroutine.  The  number  of  elements  for  the  space 
discretization,  the  length  of  the  channel  and  the  water  depth  are  first  defined.  Next  the  time  increment 
A/  is  defined  together  with  the  maximum  time  for  the  solution.  The  call  to  GEOMT  will  define  the 
particular  channel  geometry.  Following  that  step  the  user  has  to  specify  the  type  of  boundary  condition 
to  be  imposed  at  the  inlet  of  the  channel.  At  the  outlet  of  the  channel  the  boundary  condition  is 
adjusted  accordingly  by  the  code. 

NOTE:  The  number  of  elements  should  be  an  even  number.  In  the  case  of  the  half  cylinder 
geometry  this  number  should  be  a  multiple  of  twenty. 

GEOMT 

This  subroutine  is  used  to  generate  the  channel  geometry.  First  the  user  specifies  the  type  of 
geometry  desired  (more  types  can  be  added),  then  whether  the  channel  is  of  infinite  length  (open  chan¬ 
nel)  or  finite  length  (closed  channel)  is  considered.  After  the  generation  of  the  specified  geometry  a 
warning  is  given  if  the  number  of  elements  has  been  increased.  This  increase  is  necessary  if  the  spac¬ 
ing  of  the  nodes  has  been  changed  to  accommodate  a  particular  geometry.  The  different  types  of 
geometries  currently  generated  by  the  code  are  shown  in  Fig.  2. 

ELMAT 

This  subroutine  is  called  by  QSOLV  and  is  used  to  define  the  element  matrices  [A/l  and  (A],  Eq. 
(3),  and  to  assemble  the  element  matrices  [A  ]  and  (£]. 

QSOLV 

This  subroutine  performs  a  number  of  tasks.  First  the  boundary  conditions  for  the  specified  type 
are  defined.  Then  the  boundary  force  vector  is  computed  and  the  global  matrices  are  assembled  from 
the  element  matrices  defined  in  ELMAT.  A  sequence  of  calls  to  MPART,  DCOMP,  and  SOLVE  will 
partition  the  global  matrices  according  to  the  boundary  conditions  and  will  solve  the  system  of  equa¬ 
tions,  respectively. 

XOUT 

This  subroutine  is  called  by  the  main  program  and  it  is  used  for  the  output  of  results  at  selected 
devices.  The  user  has  the  choice  of  output  on  a  line  printer,  on  magnetic  tape  or  at  the  user's  terminal. 
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XPLOT 

This  subroutine  is  used  to  plot  the  results  at  selected  time  intervals  and  it  is  culled  by  XOUT. 
Results  for  the  elevation,  the  velocity  or  the  pressure  cun  be  plotted  individually  or  in  any  combination. 

TOUT 

This  subroutine  is  also  called  by  the  main  program  and  it  is  used  for  the  output  of  results  as  a 
function  of  time  and  at  specific  points  with  respect  to  the  coordinate.  Currently  the  values  of  the 
elevation  and  velocity  at  x  ■»  0  and  x  *  L  ™  (100.00)  are  obtained. 

TPLOT 

This  subroutine  is  called  by  TOUT  for  plotting  results  at  the  user's  terminal.  Both  XPLOT  and 
TPLOT  are  written  by  utilizing  the  capabilities  of  the  HP-2648  graphics  terminal. 

MPAPT 

This  subroutine  performs  the  partitioning  of  the  matrix  coefficient  of  the  system  of  equations 
according  to  the  specified  boundary  conditions. 

DCOMP  aad  SOLVE 

These  two  subroutines  perform  the  decomposition  and  solution  of  the  matrix  systems  of  equa¬ 
tions. 


The  last  three  subroutines  are  part  of  a  users  library  on  the  HP- 1000  computer  system.  This 
reduces  the  size  of  the  computer  code  and  the  compiling  time. 

4.  Pragma  Verification 

The  required  parameters  for  solving  a  test  case  with  the  developed  computer  code  are  specified  by 
the  user  through  subroutine  INPUT. 

The  example  which  will  be  used  as  a  demonstration  case  simulates  the  propagation  of  a  surface 
wave  in  a  channel  with  a  half-cylinder  placed  on  the  channel  bottom.  A  list  of  all  the  input  parameters 
as  they  appear  on  the  users  terminal  is  given  in  Appendix  A. 

The  first  parameter  is  the  number  of  elements  (NUMEL)  to  be  used  for  the  space  discretization, 
NUMEL  *  100.  The  length  of  the  computational  domain  is  XL  ~  100.00  m  and  the  water  depth  is 
Hq  «■  10.0  m.  The  length  of  each  element  is  equal  to  I  m  and  the  time  step  size  is  specified  as  Dr  * 
0.02.  The  channel  geometry  is  given  in  Fig.  3  for  the  problem  under  consideration.  For  this  type  of 
geometry  NUMEL  has  been  increased  to  1 10.  This  is  a  result  of  the  finer  mesh  where  the  half-cylinder 
is  positioned.  The  element  length  has  been  reduced  to  0.5  at  that  location. 

The  physical  problem  is  defined  by  specifying  the  type  of  boundary  conditions  at  x  —  0.0.  In  the 
present  case  the  velocity  is  specified  as 

V(0,  t )  -  V0  sin  (o»t) 

where  the  amplitude  K0  —  0.1  and  c*  ■  2 ir/T.  The  wave  period  7*  is  equal  to  5.0  which  corresponds  to 

a  wavelength  of  50  or  five  times  the  depth  Ho- 
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The  rest  of  the  program  parameters  are  needed  for  the  output  of  results  on  a  printer  and/or  tape 
and  for  the  plotting  results  on  the  users  terminal.  Examples  of  the  plotted  results  are  given  in  Figs.  4 
through  6.  In  Fig.  4  a  continuous  plot  is  given  for  the  surface  elevation.  On  the  left  side  of  the 
cylinder  one  can  observe  the  formation  of  the  standing  wave  due  to  the  wave  reflections  from  the 
cylinder.  The  actual  surface  elevation  as  well  as  the  elevation  with  respect  to  the  size  of  the  channel 
are  given  on  the  same  figure.  The  next  two  figures  are  for  the  elevation  and  velocity  at  x  -  0.0  and  x 
—  100.0.  The  change  in  elevation  due  to  wave  reflections  can  be  seen  in  Fig.  5.  These  results  have 
been  compared  to  analytical  and  experimental  ones  and  the  agreement  was  very  satisfactory  [2]. 

The  particular  problem  presented  here  is  only  an  example  of  the  type  of  problems  that  can  be 
solved  with  the  developed  code.  Additional  types  of  geometries  and  boundary  conditions  can  be  easily 

handled  by  the  code. 

IV.  CONCLUDING  REMARKS 

A  finite  element  model  has  been  implemented  for  simulating  the  motion  of  surface  waves  in  open 
channels.  The  computer  code  which  was  developed  for  that  purpose  has  been  presented  in  this  report. 
The  description  of  the  code  together  with  an  example  problem  provide  the  prospective  user  with 
sufficient  information  for  implementing  and  using  the  code.  At  this  point  the  code  can  be  used  on  any 
machine  with  a  FORTRAN  compiler.  The  only  exception  is  for  the  two  plotting  subroutines.  These 
are  designed  to  be  used  only  with  the  HP-2648  terminals.  The  code  has  been  optimized  to  a  certain 
degree  but  further  savings  in  execution  time  can  be  achieved  by  implementing  the  vector  instruction  set 
(VIS)  available  on  the  HP-1000  computer  system. 
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Fig.  1  —  Flowchart  of  computer  code 
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Fig.  3  —  Channel  geometry  for  example  problem.  Half-cylinder  with  radius  R  “  5.0  m  and  //0  -  10.0  m. 
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Fig.  4  —  Surface  elevation  along  channel  during  25  cycles  with  half-cylinder  on  the  bottom. 
Nonreflecting  boundary  at  the  right. 
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VELOCITY  5  ELEVATION 


Fig.  5  —  Time  history  of  wave  elevation  at  the  left  and  right  boundaries  for  the  case 

shown  in  Fig.  4 


TIME 

Fig.  6  —  Time  history  of  velocity  at  the  left  and  right  boundaries  for  the  case  shown  in  Fig.  4 
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****  SIMULATION  OF  WATER  WAVES  +**■* 
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****  OPEN  CHANNEL  **** 
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**************************************************  *  :*  **.*  * 
******************************************************  .** 
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****  LUIN  4  «  4,  LEFT  CASSETE  IS  THE  INPUT  DEV.  «■**» 

****  LUIN  #  -  5,  RIGHT  CASSETE  IS  THE  INPUT  DEV.  ♦  *** 

****  LUIN  #  -  8,  TAPE  DRIVE  IS  THE  INPUT  DEV.  *♦** 

****  TYPE  THE  LU  *  «  1 

*** A  **************************  ******  *******************  * 
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* *  *  *  *  *  ♦*  *  *****  **  *  **  *  *  *  *  *  *  *  *  *  *  *  *  ***  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 
****  INPUT  THE  CHANNEL  LENGTH,  < XL  »100.0>,  XL  *  100 
******************************************************** 
+***  INPUT  THE  WATER  DEPTH  HO  -  I  0 
******************************************************** 
****  INPUT  THE  TIME  STEP  SIZE  DT  -  .02 
******************************************************** 
**+*  INPUT  THE  MAX  TIME  FOR  THE  SIMULATION, TMAX  »  1. 
******************************************************** 
***«  HOW  FREQUENTLY  WOULD  YOU  LIKE  THE  RESULTS  ***<■ 

****  OUTPUTTED?  IF  IFREQ-10,  THEN  RESULTS  WILL  BE  **** 
****  OUTPUTTED  EVERY  10  TIME  STEPS.  IFPEQ  »  10 
********************************************  **  *  *  ******** 
***♦  SPECIFY  THE  TYPE  OF  CHANNEL  GEOMETRY  **** 


0043 

**** 

GEON-O 

CONSTANT  BOTTOM  SLOPE 

**** 

0044 

**** 

GEOM* 1 

CHANNEL  WITH  A  RAMP 

**** 

0045 

**** 

GEOM-2 

HALF-CYLINDER  OH  THE  BOTTOM 

*  *  ** 

0046 

**** 

GEOM-3 

BOTTOM  WITH  CYLINDRICAL  TRENCH 

**** 

0047 

0048 

**** 

**** 

GEOM  • 

GEOM-4 

2 

BOTTOM  SLOPED  UPWARD 

♦  ♦** 

0049 

***** 

******** 

I******************************* 

0050 

**** 

SPECIFY 

’  THE  TYPE  OF  CHANNEL 

**** 

0051 

**** 

TYPE-1 

INFINITE  CHANNEL 

**** 

0052 

0053 

**** 

***  + 

TYPE  - 

TYPE-2 

1 

FINITE  CHANNEL 

**** 

******************************************************** 

****  SPECIFY  THE  TYPE  OF  BOUNDARY  CONDITION  AT  INLET***’. 

***** 

****  IBC-1 ,H< 1  >«CONSTANT  **** 

****  IBC-2, H< 1  )-SINE  WAVE  **** 
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I***'?" 


0059 
0060 
0061 
0062 
0063 
0064 
0063 
0066 
0067 
0068 
0069 
0070 
0071 
0072 
0073 
0074 
0075 
0076 
0077 
0078 
0079 
0080 
0081 
0082 
0083 
0084 
0085 
0086 
0087 
0088 
0089 
0090 
0091 
0092 
0093 
0094 
0093 
0096 
0097 
0098 
0099 
01  00 
0101 
01  02 
0103 
0104 
0103 
01  06 
0107 
0108 
0109 
01  1  0 
01  1  1 
01  12 
01  13 
01  »4 
01  13 
01  16 
0117 
01  18 


****  IBC=3, V< 1  INCONSTANT  **** 

****  IBC-4, V< 1  )»SINE  WAVE  ♦*** 

****  IBC«3,H< 1 )«CONSTANT, FINITE  CHANNEL  **•* 

444*  IBC*6,H< 1  )*SIHE  WAVE, FINITE  CHANNEL  **** 

4444  IBC  ■>  2 

***%*(#**  ******  ********  ******  +  %  ************  ******  *  *  **  +  *  *  * 
4444  IF  YOU  WANT  TO  INCLUDE  THE  BOTTOM  FRICTION  *+** 

****  TERM  INTO  THE  MOMENTUM  EQUATION,  **** 

4444  TYPE  1  FOR  YES,  0  FOR  NO  :  0 

+*****i»+********»*+***++***%******+* ****** ******+******* 

****SPECIFY  THE  VALUE  OF  THE  WAVE  AMPLITUDE,  AMO  »  .1 


******************************************************** 
*#**SPECIFY  THE  VALUE  OF  THE  WAVE  PERIOD.  TO  -  3 
******************************************************** 
4444  INPUT  THE  TIME  TO  START  THE  OUTPUT  OF  RESULT  **** 
*♦**  TSTAR  ■  0. 


*4**  INPUT  THE  TIME  TO  STOP  THE  OUTPUT  OF  RESULTS  4444 
4444  TSTOP  «  I . 

TIME  «  0.000 
TIME  «  .200 

44444444444444444444444444444444444444444444444444444444 

4444  DO  YOU  WANT  TO  OUTPUT  THE  RESULTS  ON  THE  444* 

*444  PRINTER?  TYPE  1  FOR  YES,  0  FOR  NO  :  0 

4444444444444444444444444444444444444444*444444444444444 

4444  DO  YOU  WANT  TO  OUTPUT  THE  RESULTS  ON  THE  TAPE  **** 
4444  DRIVE?  TYPE  1  FOR  YES,  0  FOR  NO  i  0 


4444SPECIFY  WHICH  OF  THE  FOLLOWING  VARIABLES  YOU  ***4 
4444WANT  PLOTTED  i  4444 
4444  0  m  NOTHING  4444 
*444  1  «  ELEVATION  **** 
444*  2  «  VELOCITY  4444 
444*  3  a  PRESSURE  4444 
4444  4  -  ELEVATION  AND  VELOCITY  **44 
4444  3  a  VELOCITY  AND  PRESSURE  *•*** 
4444  6  4  ELEVATION  AND  PRESSURE 

*444  7  ■  ELEVATION. VELOCITY, AND  PRESSURE  "*•* 
4444  I PLOT  ■  0 


TIME  ■  .400 
TIME  «  ,600 
TIME  ■  ,800 
TIME  •  I . 000 

NUMBER  OF  POINTS,  N P  •  111 


THE  «  OF  POINTS  FOR  THE  TIME  PLOT  IS  -  3 

*4*444444444444444444444444  444444444  *4  44  4 4*4 **4444 444444 

4444  DO  YOU  WANT  PRINTER  OUTPUT  OF  THE  NUMBER  OF  4444 
4444  POINTS, PLOTS, ETC.?  TYPE  1  FOP  YES,  0  FOR  NO  :  0 
*4*4444444444444444444444444444444444444**444444**4*4*4* 
*444  00  YOU  WANT  TAPE  OUTPUT  OF  ELEVhTIOM  AND  *♦♦* 

4444  VELOCITY?  TYPE  1  FOR  YES,  0  FOR  NO  0 
READY  TO  PLOT  GRAPHICS,  PPESS  RETURN 
/ 

**4*444*4*4444*4  44**  4  *•*  *  **  *4444*  ***  4  *  4  **  -  -  ,•*«*********-, 


ii 


0001  FTN7X, L 


0002  PROGRAM  KYMA<3,93),  UPDATED  9:20:92 

0003  £***** 

0004  C*********************— *************************************** 

0005  C*****  ***** 


0006  C*****THIS  PROGRAM  SOLVES  THE  DEPTH  AVERAGED  CONTINUITY  ***** 

0007  c*****AND  MOMENTUM  EQUATIONS  IN  ONE  DIMENSION.  THE  ***** 

0009  C*****C9MPUTATI0NAL  MODEL  IMPLEMENTED  IN  THIS  PROGRAM  ***** 

0009  C*****IS  BASED  ON  THE  FINITE  ELEMENT  METHOD.  THE  MODEL  ***** 

00(0  C*****IS  DERIVED  BY  ASSUMING  LINEAR  APPROXIMATION  IN  ***** 

0011  C*****SPACE  AND  TIME.  ***** 

0012  C*****  FLUID  DYNAMICS  BRANCH  ***** 

0013  C*****  NAVAL  RESEARCH  LABORATORY  ***** 

0014  C*****  ***** 

0015  C* ************************ ****** ******************************* 
0016  C****» 

0017  COMMON/BLOC 1 /  TO, T, TMAX, I FREQ, IBC, NUMBC, IB<  2  > 

0018  C0MM0N/BL0C2/  G< 222),  YQ< 222 > 

0019  C0MM0N/BL0C3/’  H<  1 20  ),  V<  1 20  >,  DH<  1 20  >,  DV<  1 20  >,  PR<  1 20  > 

0020  C0MM0N/BL0C4/  RH< 1 20 > , RL< 1 20 >, DX< 1 20 > , HC< 1 20 > 

0021  C0MH0N/BL0C5/  N,WO,DT,GE,HO,XL,CZ,  IGEOM,  ICZ.AMO 

0022  C0MM0N/8L0C6/  HTN< 151 ), VTN< 151 >,HT1< 151 >, VT1< 151 >, TIN< 151 ) 

0023  DIMENSION  GSM< 2503), IP<222> 

0024  C***** 


0023 

WR I TE< 1,1000) 

0026 

*** 

0027 

1  ,/,-**** 

****" 

0028 

2  ,/,«**— 

0029 

3  ,/,«**** 

COMPUTATIONAL  HYDRODYNAMICS 

***** 

0030 

4  ,/,"**** 

***** 

0031 

5  ,/,•**** 

SIMULATION  OF  WATER  WAVES 

***** 

0032 

6  ,  /',■**— 

IN  AN 

***** 

0033 

7  ,/',-**** 

OPEN  CHANNEL 

***** 

0034 

8  ,/,"***• 

***** 

0033 

9  . / , ■ ******* 

*** 

0036 

GO  TO  5 

0037 

100  CONTINUE 

* 

0038 

WRITE< 1,1040) 

0039 

WRITE< 1.1003) 

0040 

1003  FORMAT< ******** 

*** 

0041 

1  *****  00 

YOU  WANT  TO  RUN  ANOTHER  CASE? 

***** 

0042 

2  ,/,-****  TYPE 

1  FOR  YES,  0  FOR  NO  :  ) 

0043 

READ< 1 , *  >  ICASE 

0044 

IF< ICA8E.LT. 1 > 

GO 

TO  9999 

0043 

WRITE< 1,1010) 

0046 

1010  FORMAT< ■*•**•** 

+++ 

0047 

1  ,/,-****  WOULD 

YOU  LIKE  TO  KEEP  SOME  OF  THE  DATA  THE 

***** 

0048 

2  ,✓,«****  SAME? 

TYPE  1  FOR  YES,  0  FOR  NO  :  _“ > 

0049 

READ<  LUIN, * )  IQUES 

0050 

IF< IQUES.LT. 1 ) 

I CALL- 0 

0051 

IF< IQUES.LT. 1 ) 

GO 

TO  5 

0052 

WRITE< 1,1013) 

0053 

*•* 

0054 

1  ,/,-****  WHICH 

SET  OF  DATA  WOULD  YOU  LIKE  TO  CHANCE? 

***** 

0055  2  ,✓,■-***  O-MORE  THAN  ONE  SET  ****• 
0056  3  ,✓,»****  1-LU  *  FOR  INPUT  OF  DATA  ****" 
0057  4  ,/,-****  2-INPUT  DATA< CHAHNEL  LENGTH, TIME  STEP. ETC.  )  ♦**♦" 
0038  5  ,/,«****  3-CHANNEL  GEOMETRY  ***** 


12 


6  ,/,•****  WHICH  SET?  ") 

READ< LUIN, - >  I CALL 
IF< ICALL .  GT .  1 >  GO  TO  6 
9  WRITE< 1 , 1 020 > 

t  020  FORMAT< ■************************************************, 

1  ,/,*****************  LU  FOR  INPUT  OF  DATA  **********. 

2  ,✓,■****  LU1N  •  ■  1,  TERM  IS  THE  INPUT  DEVICE 

3  ,/,■****  LUIN  •  -  4,  LEFT  CASSETE  IS  THE  INPUT  DEV. 

4  ,/,”*•**  LUIN  •  ■  9,  RIGHT  CASSETE  IS  THE  INPUT  DEV. 

3  ,/,«****  LUIN  *  «  8,  TAPE  DRIVE  IS  THE  INPUT  DEV. 

6  ,/,■****  TYPE  THE  LU  *  »  _“ > 

READ<  t ,*>  LUIN 

IF< ICALL. EQ. I >  GO  TO  7 


***** 

***** 

***** 

***** 


CALL  INPUT< LUIN, ICALL) 
IFCN.EQ. 0)  GO  TO  9999 


PI-3. 14199 
GE-10.0 

C0-SQRT<  HO-GE > 

NUHEL-N 

NPOIN-N+1 

NEQ-2-NP0IN 


I GRAVITATIONAL  ACCELERATION 
•WAVE  VELOCITY 

• NUMBER  OF  ELEMENTS 
•NUMBER  OF  NODAL  POINTS 
•NUMBER  OF  EQUATIONS 


C***«* INITIALIZE  PROGRAM  PARAMETERS 
C— — 

T-0.0 

DO  10  1-1, NEQ 

IP<I)-0.0  I  VECTOR  FOR 

Q< I )-0. 0  'VECTOR  OF 

YQ< I )-0 . 0  'VECTOR  OF 

10  CONTINUE 

DO  20  1-1 , HPOIN 
V<I)-0.0  I  VELOCITY  V 

H< 1  )— 0 . 0  (ELEVATION 

RH<  I  )-HC<  I  )  !  TOTAL  WATEI 

PR< I )-0 . 0  (PRESSURE  V 

20  CONTINUE 


I  VECTOR  FOR  PI VOTING< DCOMP, SOLVE > 
•VECTOR  OF  NODAL  UHKNOUNS< TIME  T> 
•VECTOR  OF  NODAL  UNKNOWNS< TIME  T-DT ) 


(VELOCITY  VECTOR 
(ELEVATION  VECTOR 
(TOTAL  WATER  DEPTH  VECTOR 
•PRESSURE  VECTOR 


20 

C*«*** 


NBAND-7  (BANDWIDTH  0 

NHBAN-<  NBAND- 1  V2 
JGF-NHBAN-NEQ 
JGSM- JCF ♦NEQ 
JEND-JGSM+NEQ-NBAND 
IF< JEND.CT.2909>  CO  TO  8888 
1CUNTE-0 
IT-0 
IPT-0 
HPLOT-O 


BANDWIDTH  OF  MATRIX  EQUATION 


WRITE< 1,1029) 

1  029  FORMAT<  ■••********••***•**•*•— ************************** 

1  ,/,**••*  INPUT  THE  TIME  TO  START  THE  OUTPUT  OF  RESULT 

2  ,/,•****  TSTAR  -  ") 

READ< LUIN, - )  TSTAR 
WR1TE< 1,1030) 

1  030  FORMAT<  •****-*********—****♦*—*****♦***  ****♦***♦*.«.****« 
t  ,/,-****  INPUT  THE  TIME  TO  STOP  THE  OUTPUT  OF  RESULTS 


***** 


0119  2  TSTOP  -  " > 

0120  READ< LUIN, * >  TSTOP 

0t2t  190  CONTINUE 

0122  WRITE< 1 , 1 039  >  T 

0123  IFCT.LE.TSTAR. OR. T.GE. TSTOP)  CO  TO  200 

0124 

0125  CPU.  XOUT< LUIN, IPT  > 

0126  C***** 

0127  C***1* '•STONE  N<  0,  T),H<  L,  T  ),  V<  0,  T>.  V<L,  T  ) 

0120  C***** 

0129  IT-IT-M 

0130  NTN<  XT )«H< NPOIN >  IM<L,T> 

0131  VTN< IT >-V< NPOIN)  »V<L,T> 

0132  NT1<  IT)«H<  1  )  IH<0,T> 

0133  VT1<  IT )«V<  1  >  IV<0,T) 

0134  TIH< IT  >«T 

0139  IF< IT.CE. 151 >  GO  TO  46 

0134 

0137  200  T-T+OT 

0139  C***** 

0139  C*****TNE  ARRAY  C9N< I )  IS  USED  FOR  STORING  THE  N0AND 

0140  C**»**DIAGONALS  OF  THE  ST1FNESS  MATRIX 

0141  C*«***GSH< I >.  I-JGF+1 , JG9M,  RIGHT  HAND  SIDE  VECTOR  OR 


0142  O****  SOLUTION  VECTOR. 

0143  C«*«**GSIKI),  I-JCSN-H  ,  JEND,  1  ST ,  2ND ,  ,  7TH  NON- 

0144  ZERO  GIAGONALS . 

0145 

0146  DO  30  1-I.4END 

0147  30  G8N<  I  )«0. 0 

0148 

0149  CALL  QSOLV<  GSM  I  >,  VEND,  I P , NEQ , NDAND ) 

0150  C***** 

0151  DO  40  1-1. NCQ 

0152  Q< I )"CSM< I ♦ JGF ) 

0153  40  CONTINUE 


0154  C***** 

0199  C«**«*CONPUTE  THE  ELEVATION.  VELOCITY  4  PRESSURE 
0156  C***** 

0157  DO  45  I«1. NPOIN 

0158  H<I)«0<2*I-1  ) 

0159  V< 1 )*G< 2* l ) 

0160  RN< I >«MC< I >*M< I ) 

0141  PR<!>-H<!)*0.S*V<I)**2 

0162  45  CONTINUE 

0163  C**+** 

0164  C*****IF  T  EXCEEDS  THAN,  TERMINATE  INTEGRATION.. 

0145  €*•**«• 

0146  IF<  T . LE. TMAX)  GO  TO  50 

0147  €*•<•«** 

0148  CO  TO  49 

0169  44  UR1TE< 1.47) 

0170  47  FORMAT<  •*•*•*••••**••••*****  WARNING  •<****■ 


0171  1  ,✓,•*•**  THE  NUMBER  OF  POIHTS  FOR  THE  TIME  PLOTS  HAS  *♦■**" 

0172  2  ./."••**  EXCEEDED  THE  MAXIMUM  NUMBER  OF  POINTS  ALLOWED.  •+*** 

0173  3  ./,-**•*  INCREASE  IFREQ  FOR  LESS  TIME  PLOTS  IF  ANOTHER 

0174  4  ,/,••***  CASE  IS  DESIRED.  ***<•" 

0175  5  ./. -*•••*•••••••••**+•••***•*+**•**•**•••**•••••**••*•*••*••• > 

0176 

0177  49  CALL  TOUT<LUIN. NPOIN, IT, IPT > 

0178  C*4*** 
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0179 
0190 
otat 
0182 
0193 
0194 
0199 
0186 
0187 
Ot  98 
0189 
»  0190 

0191 
0192 
0193 
0194 
01 93 
0196 
0197 
0198 


4KYMA 

0199 
0200 
0201 
0202 
0203 
0204 
0203 
0206 
0207 
0209 
0209 
0210 
021 1 
0212 
0213 
0214 
0219 
0216 
0217 
0219 
0219 
9220 
0221 
0222 
0223 
0224 
0223 
0226 
0227 
0228 
0229 
0230 
0231 
9232 
0233 
0234 
9239 
0236 
0237 


GO  TO  100 

50  ICUNTE*ICUNTE+1 
00  33  I>1 , NEQ 
YQ<  I  )>Q<  1  > 

33  CONTINUE 
C»**> 

C . . PRINT  RESULTS  OR  CONTINUE.. 

C*»** 

IF  < ICUNTE.NE.  I FREQ )  GO  TO  200 
1CUNTE-0 

GO  TO  130 
8888  CONTINUE 
9999  STOP 

1033  FORMAT  <//',3X,6HTIME  -,F6.3,//> 
104  0  FORMAT*  > 

END 


T-00003  IS  ON  CR00011  USING  00168  BLKS  R-0000 


c************************************************ 

C*M*» 

SUBROUTINE  INPUT<  LUIN, ICALL  >,  N.C.CHU 

C*»*> 

COMMON/BLOC 1  /  T  0 ,  T , TNAX . I FREQ , I BC , NUM8C , IB< 2  > 

C0MM0N/BL0C2/  Q< 222 > , YQ< 222 > 

COMMON/BLOC3/  H<  1 2 0  ) .  V<  1 2 0  > ,  DH<  1 2 0  > ,  DVO  2  0  > ,  PR<  1 2  0  > 

C0MM0N/8L0C4/  RM< 1 20  >, RL< 1 20  > . OX< 1 20  > . HC<  1 20 ) 

C0HM0N/BL0C3/  N, WO.OT.GE.HO.XL.CZ,  IGEON,  ICZ.AMO 

C***« 

IF< ICALL.LT. 1 )  GO  TO  10 
WRITE* 1 , 1 000 )N 

1 000  FORMAT* ••••*******i|m«**«*im*************»*»**.i******************»" 

1  ,/,■••••  HIT  THE  CARRIAGE  RETURN  FOR  ANY  VARIABLE  THAT 

2  ,/,«***•  YOU  WOULD  LIKE  TO  REMAIN  THE  SAME.  »»• 

3  ■**•*>>>>>>>>>>>>>>>>>  WARNING  <<<<<<<<<<<<<<<<<<<<»»• 

3  ,✓,•♦***  THE  ELEMENT  NUMBER  HAS  BEEN  CHANGED  TO  “,I4,"  **»“  > 

10  WR1TE< 1 . 1010) 

1010  FORHAT< ■•••***•*•*•*******#*******•*****************************" 

1  ,/,"**•*  INPUT  NUMBER  OF  ELEMENTS.  IF  A  HALF  CYLINDER  •»** 

2  ,/,•****  IS  GOING  TO  BE  USED  FOR  THE  CHANNEL  GEOMETRY.  ****• 

3  ,/.-*•**  THE  NUMBER  OF  ELEMENTS  SHOULD  BE  A  MULTIPLE  OF  ****- 


4  •*»**  20.  NELEM  •  "  > 

READ*LUIN,*>  N 

IF< ICALL. EQ. 3)  GO  TO  20 

WRITE< 1.1013) 

1 013  FORMAT< ■••*#•••♦••**♦•*•***••****♦•***♦*******•*< .**•**•***»**•**" 

1  THE  LENGTH  OF  EACH  ELEMENT  L  IS  DEFINED  AS  #*•*- 

2  DX<  L  >  ■  XL/NEMEL  »»« 

3  WHERE  XL  IS  THE  LENGTH  OF  THE  CHANNEL 

4  , /.  ■*•••••**>•*•*••****#****♦*»*****.*****»*»*.*****♦♦>**»»  «  > 
HRITE< 1,1020) 


1  020  FORMAT<  ■•••*•••••••••••♦••**••••**•*•*••••*.*•*****••*•*•••****••* 

I  ,/,>**•*  INPUT  THE  CHANNEL  LENGTH,  <XL  >100.0),  XL  •  * > 

RtAO<LU!N,«>  XL 
20  WO«PLOAT<  N  VXL 
HPOIH>H+1 

IF< ICALL. EO. 3)  GO  TO  30 
WRITE<  1 , 1023) 
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0238  1025  FORNAT<  ********************************************************** 

0235  1  ,/.•****  INPUT  THE  WATER  DEPTH  HO  -  _“  > 

0240  READ< LUIN,*)  HO 

0241  WRITE< 1 , 1 030 > 

0242  1 030  FORNAT<  ********************************************************** 

0243  1  ,✓,•****  INPUT  THE  TIME  STEP  SI2E  DT  -  > 

0244  READ<  LUIN, *>DT 

0245  WRITE< 1.1035) 

0246  1 035  FORMAT< ********************************************************** 

0247  1  ,✓,•****  INPUT  THE  MAX  TIME  FOR  THE  SIMULATION. TMAX  -  - > 

0246  READ<LU1N,«)  TMAX 

0249  WRITEO,  1040) 

0250  1  040  FORHAT<  ********************************************************** 

0251  1  ,/,•****  HOW  FREQUENTLY  WOULD  YOU  LIKE  THE  RESULTS  ***** 

0252  2  ,/,•****  OUTPUTTED?  IF  I FREQ-10,  THEN  RESULTS  WILL  BE  ***** 

0253  3  ,/,-****  OUTPUTTED  EVERY  10  TIME  STEPS.  I FREQ  »  _• > 

0254  REAO< LUIN,*>  I FREQ 

0255  C***«* 

0256  C«****SET  CHANNEL  GEOMETRY 


0257 

0258 

0259 

0260 

0261 

0262 

0263 

0264 

0265 

0266 

0267 

0268 

0269 

0270 

0271 

0272 

0273 

0274 

0275 

0276 

0277 

0276 

0279 

0280 


0263 

0264 

0265 

0266 

0267 

0288 

0289 

0290 

0291 

0292 

0293 

0294 

0295 

0296 

0297 

0296 


30  CALL  GEOHT< NPOIN.LUIN, IT, ICALL) 

C***** 

IF< ICALL • EQ . 3 )  GO  TO  40 
WRITE< 1,1045) 

1 045  FORHAT< •********************************************************“ 
1  ,/,•****  SPECIFY  THE  TYPE  OF  BOUNDARY  CONDITION  AT  INLET***** 


2  ,/,“*•**  ***** 

3  ,/,•****  IBC-1 ,H< 1 )-COHSTANT  ***** 

4  ,/,**••*  IBC-2.W 1 )«SINE  WAVE  ****« 

3  ,/,«**•*  IBC-3, V< 1 >-COMSTANT  ***** 

6  ,/,•****  IBC-4, V< 1 )-SINC  WAVE  ***** 

7  ,/,-*•**  1BC-S,H< 1 )-C0N8T ANT, FINITE  CHANNEL  ***** 

8  ,/,-*•••  IBC-6,H< 1 )-SINE  WAVE, FINITE  CHANNEL  ***** 

9  ,/,■**•*  IBC  -  _*> 

REAO< LUIN, • )  IBC 

40  CONTINUE 


NUH6C-2  !  NUMBER  OF  BOUNDARY  CONDITIONS  REQUIRED 

GO  TO  <30,50,51,51,50,50),  IBC 

50  IB< 1 >*1  I  INDEX  FOR  BOUNDARY  CONDITION  AT  X-0 

IB<  2  >-2*<  N*1 )  »  INDEX  FOR  BOUNDARY  CONDITION  AT  X-100 
GO  TO  57 

51  IB< 1 >*2 
16<2)-2*<N+I  )-1 

57  1F< I CALL. EQ. 3)  GO  TO  60 
WRITE< 1 , 1 050) 

1 050  FORHAT<  ********************************************************** 

1  ,/,•****  IF  YOU  WANT  TO  INCLUDE  THE  BOTTOM  FRICTION  ***** 

2  ,/,•****  TERM  INTO  THE  MOMENTUM  EQUATION,  ***** 

3  ,✓,■*•**  TYPE  1  FOR  YES,  0  FOR  NO  i  _"  > 

READCLUIN,*)  ICZ 

IF< ICZ.LE . 0)  GO  TO  59 
WR I TE< 1,1055) 

1 055  FORMAT< ********************************************************** 
I  ,/, *****SPECIFY  THE  VALUE  OF  THE  CHEZY  COEFFICIENT.  CZ  •  > 

READ< LUIN, * )  CZ 
59  WRITE< 1 , 1 060  ) 

1 060  F0RMAT< »********************************************************" 
1  *****8PECIFY  THE  VALUE  OF  THE  WAVE  AMPLITUDE.  AMO  •  > 

Rf AO< LUIN, * )  AMO 
CO  TO  <60,70,80,70,80,70), IBC 
70  WRITE< 1,1065) 
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0299  1 065  FORMAT< •----4*--4*-4--****-****-*4-*-*4**-4>--***-****-**** 

0300  1  ,/, “444-SPECIFY  THE  VALUE  OF  THE  WAVE  PERIOD.  TO  ■  "> 

0301  READ<  LUIH, * )  TO 

0302  90  RETURN 

0303  END 

0304 


4.KVNA  T-00003  IS  ON  CR00O1 1  USING  00169  BLKS  P-0000 

0305  C*’*****4*4***4— •**—«*—*— *♦—♦-♦****-♦*— *  — 

0306  C— — 

0307  BLOCK  DATA 

0308  COHHON/BLOC1/  TO, T, THAX, I FREQ, IBC, NUMBC, IBC  2  > 

0309  C 0HH0N/9L 002/  Q< 222 >, Y0< 222 > 

0310  C0HH0H/BL0C3F  H<  1 20  > ,  V<  1 20  > ,  DH<  1  20  > ,  DV<  1 20  > ,  PPM  20  > 

0311  C0NH0N/BL0C4/  RN< 1 20 >, RL< 1 20 >, 0X< 1 20 > , MC< 1 20 > 

0312  COriNON/BLOCS/’  N,W0,DT,GE,H0,XL,CZ,  IGEOM,  ICZ,hMO 

0313  C0NN0N/BL0C6/  HTN<  131  >,VTN<  131  >,HT1<  131  >,VTK  131  >,TIM<.  131  > 

0314  END 

0313  C«—*  * 
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0316  c**«*******4**<*-— ********* ************ ******** 

0317  C***** 

-0318  SUBROUTINE  GEONT< NPOIN,LUIN, IT, ICALL),  N.C.CHU,  G. A.KERAMIDAS 

0319  C***** 

0320  C0NH0H/BL0C3/  M<  1 20  >,  V<  1 20  ), OH<  1 20  >,  DV<  1 20  >,  PR<  1 20  » 

0321  COHHONFBLOC4/'  RN< 120), RL< 120 ), DX< 1 20  >, HC< 1 20  > 

0322  COHNON/BLOCS/  N.WO,DT,CE.H0,XL.CZ, ICEOH, ICZ,ANO 

0323  C***** 

0324  IF< ICALL. EQ .2)  GO  TO  3 

0323  WRITE< 1,1 000  > 

0326  1  000  FORMAT*  'mm****************************************************' 

0327  1  ,/,•****  SPECIFY  THE  TYPE  OF  CHANNEL  GEOMETRY  ***•" 

0328  2  ,/,•****  GEON-O  CONSTANT  BOTTOM  SLOPE  —  **• 

0329  3  ,/,«•***  GEOM-1  CHANNEL  WITH  A  RAMP  **** ■ 

0330  4  ,/,•****  GEOM-2  HALF-CYLINDER  ON  THE  BOTTOM  ***•■ 

0331  3  ,/,•****  GEOM-3  BOTTOM  WITH  CYLINDRICAL  TRENCH  •***■ 

0332  6  ,/,•****  CEOM-4  BOTTOM  SLOPED  UPWARD  —  *• 

0333  7  ,/,-**•*  CEON  -  ■ > 

0334  READ<LUIN,*>  IGEOH 

9333  WRITE< 1 , 1 01 0> 

0336  1010  FORHAT<  ********************************************************** 

0337  1  .A"****  SPECIFY  THE  TYPE  OF  CHANNEL  ****  • 

0339  2  ,/,'••••  TYPE-1  INFINITE  CHANNEL  ****  “ 

0339  3  ,/,■*•••  TYPE-2  FINITE  CHANNEL  ***♦" 

0340  4  ,/,•****  TYPE  -  > 

0341  READ< LUIH, * >  I TYPE 

0342  IF< ITYPE.E0.1 >  GO  TO  3 

0343  WR1TE< 1,1013) 

0344  1 013  FORNAT< •***•••****•*•****••*** **************************** ******* 

0343  1  ,/,•****  INPUT  THE  DEPTH  OF  THE  CHANNEL  AT  THE  CLOSED  **•*• 

0346  2  ,/,•****  END,  HNO  ■  " > 

0347  READ< LUIH,*)  HNO 

0349  GO  TO  6 

9349  3  CONTINUE 

0330  IP< ICEON, HE, 4 >  GO  TO  6 

0331  WRITK  1 ,  1920> 


0352 

1020 

FORHAT< “ ********************************* *************** ******** 

0353 

1 

,/,-*•*•  INPUT  THE  DEPTH  OF  THE  WATER  AT  THE  END  OF  **** 

0354 

2 

;  ,/,«****  THE  CHANNEL  DWO  ■  "> 

0355 

R£AD<  LUIH, * )  DUO 

0350 

HNO-HO-OWO 

0357 

6 

DXO*t . O/WO  »  CONSTANT  ELEMENT  LENGTH  *  XL/N 

0358 

DO  7  1*1 , NPOIN 

0359 

DX< I >-0X0  !  ELEMENT  LENGTH  VECTOR 

0360 

HC< I  )*H0  !  WATER  DEPTH  IN  CHANNEL 

0361 

RL< : >»DX< I >*FLOAT< 1-1  )  *  NODAL  POINTS  COORDINATES 

0362 

7 

CONTINUE 

0363 

IF< I TYPE . EO . 1 >  GO  TO  9 

0364 

HC<NP01N>— HNO 

0365 

9 

IF< IGEOH.LE. 0>  GO  TO  50 

0366 

GO  TO  < 1 1.21,31 ,39  >,  IGEOM 

0367 

c***—* 

0368 

c> 

'GEOMETRY  FOR  CHANNEL  WITH  A  RAMP 

0369 

c**-** 

0370 

1 1 

CONTINUE 

0371 

HN  0*5 

0372 

NEH-N/2+1 

0373 

NEC-N/1 0 

0374 

NED-N£H*NEC 

0375 

DO  15  I-HEH.HED 

0376 

DX< I >-0X0/2 . 0 

0377 

15 

HC< 1 )— HO— FLOAT< I-NEH >*DX< I > 

0378 

NMAX-HPOIN+NEC 

0379 

DO  20  I-HED.NMAX 

0380 

DX< I >-0X0 

0381 

20 

HC< I >— HO-HHO 

0382 

GO  TO  43 

0383 

c< 

1 

0384 

c*****geometry  for  channel  with  a  half-cylinder  on  the  bottom 

0385 

c****« 

1  - 

0386 

21 

HCC-H/10 

0387 

NCO-<  N— NEC  V2 

0388 

NCI -HC 0*1 

0389 

NC2-NC  0*2*NEC 

0390 

DR— XL/1 0 . 0 

0391 

RR-DR/2 . 0 

0392 

C’ 

0393 

DO  25  1-1, NCO 

0394 

DX<  I  >— DXO 

0395 

DX<  I*NC2  >-0X0 

0396 

25 

CONTINUE 

0397 

c< 

•398 

DO  27  1—1 ,NPOIN*HEC 

•399 

NC<  I  >— HO 

0400 

27 

CONTINUE 

0401 

C< 

1 

0402 

DO  30  I-NC1.NC2 

0403 

11-1-<NEC/2*1 > 

•404 

DX< I >—0X0/2 . 0 

•405 

XR— XL/4 . 0*DX<  I  >*FLOAT<  1 1  > 

0406 

HC< 1 >-N0-08RT<  RR**2-XR**2 > 

•407 

30 

CONTINUE 

0408 

c 

***** 

•409 

CO  TO  43 

0410 

C 

***** 

1 

04 11 

C*****CEOHETRV  FOR  A  CHANNEL  WITH  A  CYLINDRICAL  TRENCH  ON  THE  BOTTOM 

0412 

C< 

***** 

1 
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0413 

31 

HEC-N/10 

0414 

NC0«<  N-NECV2 

0419 

NC1-NC0+I 

•416 

HC2-NC 0+2*N£C 

0417 

DR-XL/10.0 

0419 

RR-DR/2. 0 

0419 

C* 

» 

0420 

00  36  X-1,NC0 

0421 

DX< I  >-DX0 

•422 

OX(  I*NC2>«DX0 

0423 

36 

CONTINUE 

•424 

O 

» 

0425 

DO  37  1-1 .NPOIN+NEC 

0426 

HC< I >-N0 

0427 

37 

CONTINUE 

0429 

O 

***« 

1 

0429 

DO  39  I-NC1.NC2 

0430 

I 1-I-< NEC/2+1 > 

0431 

OX< I >-0X0/2 . 0 

0432 

XR— XL/4 . 0*OX<  I  >-FLOAT<  1 1  > 

0433 

HC<  I  >-H0*S0RT<  RR—2-XR—2  > 

0434 

39 

CONTINUE 

0433 

GO  TO  43 

0436 

c***** 

•437 

C«**— GEOMETRY  FOR  A  CHANNEL  WITH  THE  90TT0M  SLOPED 

UPWARD 

0439 

C* 

» 

0439 

39 

SLOPE -HN9/XL 

0440 

DO  41  1—1 , HPOIH 

0441 

HC< I >-N0— OX< I >*SL0PE-FL0AT< 1-1 > 

0442 

41 

CONTINUE 

0443 

c— — 

0444 

43 

M-M+MCC  1  THE  MUH9ER  OF  ELEMENTS  HAS 

9EEH  INCREASED 

0443 

1 

C  !  5Y  NEC  FOR  THE  SPECIFIED  GEOMETRY 

0446 

DO  43  1-1, N 

0447 

RL<  1*1  >-RL<  1  >+DX<  I  > 

•449 

43 

CONTI HUE 

0449 

C* 

• 

•430 

30 

CONTINUE 

•451 

RETURN 

•452 

EHO 

•453 
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•454  C****— 1 — — ■ 

•455  C*»»* 

•454  SUBROUTINE  ELHAT < L , ESN  1 , CtN2 > .  G. A.KERAHIOAS 

•457  C*m** 

•455  C*—— SUBROUTINE  ELHAT  EVALUATES  TNI  KLCNCNT  MATRICES 

•459  f*»— 

9449  C0HH0N/9L0C1/  TO, T,THAX,  XFRCQ,  X9C.NUH0C,  I5<  2  > 

•441  COMMON/ 9L0C 3/  N<  12«>,V<  120>.DH<  12»>,DV<  I2»>,PR<  120> 

9443  C0HH0N/5L0C4/  RM<  1 24 >, RL<  120 >, DX<  1 20  ),HC<  (20  > 

•4«3  COHHOM/9LOC8/  N,tl«.0T.CC,H*,XL.C2, IG10H, ICZ.AHO 

•444  OINCNtlON  M  4, 4  >.  •<  4,  4  >,  CSH1<  4, 4  >,  BSH2<  4, 4  ) 

•445  (****• 

•444  C****MlM9t  N4TKIK,  CORRESPONDS  TO  MATRIX  CHJ  OF  CO  ■  <3> 
•447  €••••* 

•445  A<t,1>«2./DT 

M l,3>-1 ./or 

4<  ».4>«0.0 


•449 

9479 

•471 
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0472 

0473 

0474 

0475 

0476 

0477 

0475 

0475 

0450 

0451 

0452 

0453 

0454 

0455 

0456 

0457 

0485 

0488 

0480 

0481 

0482 

0483 

0484 

0485 

0486 

0487 

0498 

0499 

0500 

0501 

0502 

0903 

0504 

0505 

0506 

0507 

0508 

0509 

0510 

0511 

0512 

0513 

0514 

0515 

0516 

0517 

0518 

0519 

0520 

0521 

0522 

5523 

0524 

0529 

0526 

0527 

0528 

0529 

0530 

0531 

0832 

0933 

0834 


A<2, I >>0.0 
A<2.2»2./0T 
A<2,3>>0.0 
A<2.4>>1  •  SOT 
A<  3. 1  >>1  ./OT 
A<3,2>-0.0 
A<3.3»2./DT 
A<  3.  4  >>0 . 0 
A<  4, 1 >>0.0 
A<  4. 2  >>1 ./or 
A< 4*3 >>0.0 
A<4*4>-2./DT 

C»*» 

c**»*convective  matrix,  corresponds  to  MATRIX  CKJ  OF  EQ.  < 3 > 
c>»» 

C0>S8RT<  GE*RH<  L  >  > 

DV0>V<  L  )— V<  1  > 

DIF  0>SQRT<  2 . *OX<  L  > A  C  0>DT  >  > 

DIF1 >1 ,/OIFO 
OXF2>OIFO>*2 

DIFL><  DIF2*AB8< DVO  )+C 0*OIF 1 >/DX< L > 

C*»» 

HL  0-0 . S*<  HC< L >+HC< L* 1 >  > 

HL 1  -<  2 .  *HC<  L  >>HC<  L*  1  >  VDX<  L  > 

HL2«<  HC<  L  >«-2  .  >HC<  L ♦  1  >  VDX<  L  > 

DHDX><  M<  L*  1  >-H<  L  >  >/DX<  L  > 

DVDX><  V<  L+ 1  >-V<  L  >  VOX<  L  > 

OHC><  HC<  L*1  >-HC<  L  >  VOX<  L  > 

AVEL-0 . 5><  V<  L  >*V<  L*  1  >  > 

AOEP-0 . 5*< RH< L  >+RH<  L*t  >  > 

D ISP-0.0 

IF<  IC2.GT.  0>  DISP-GE/<  ADEP*CZ*'*2  > 

D1S1 1 -018FX  V<  L  >«AV1L  > 

DIS22-DISPX  V<L*1  >+AVEL> 

DI312-0ISP*AVEL 

DIS21-DISP*AVEL 

C***» 

8<  2.2>>2.>DVDX*DIFL+DIS1 1 
B<  2 « 1  »3 .  >G£/DX<  L  > 

5<2.4»OVOX-OIFL*DISt2 

8<2.3»3.*GE/0X<L> 

8< 1  *  2  >>2 . *<  DMOX+OMC  >-ML 1 

8<  1 , 1  »2 .  *DVDX 

8<  1 .4»0HDX*DHC*HL1 

5<I,3»0VDX 

8<4.2>*©VDX-DIFL*DIS21 

K4,l  >— 3.*CE/DX<L> 

5<  4*  4 >»2 . *0VDX*DIFL+DIS22 
5<4,3>— 3.*CE/0X<L> 

*3,2 >-0HOX+DHC-HL2 
8<3.1 >-OVOX 

8<  3 . 4  »2 .  X  OHOX+DHC  »»ML2 
K3.I  >>2  .•DVDX 

DO  10  I>! .4 
DO  10  >1.4 

o**>* 

C*»»8LtMtNT  STIFFNESS  MATRICES. 

C«>**> 

ESM 1 < I , J >>A< I ,  J  >♦  0 . 5*B< I , J  )  ’MATRIX  CAJ  OF  EQ . < 4 > 

ESM2< I , J »A< I , J >-0 . 5>8< I . J  >  ! MATRIX  C81  OF  EQ . ( 4  ) 

10  CONTINUE 
C**4M 


20 


A 


0933  RETURN 

0936  EHO 

0337  O'**** 
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0938 

0339 

0940 

0341 

0942 

0943 

0344 

0949 

0346 

0347 

0949 

0949 

0330 

0991 

0992 

0993 

0394 

0933 

0996 

0997 

0998 

0399 

0960 

0961 

0362 

0963 

0964 

0969 

0966 

0367 

0968 

0369 

8970 

0971 

0972 

0973 

0974 

0973 

0976 

0377 

0978 

0979 

0990 

0991 

0992 

0993 

0994 

0999 

0396 

0997 

0399 

0399 

0390 

0991 

0992 

0993 

0994 

*395 


SUBROUTINE  GS0LV< GSM, 4END, IP, NEQ, NBAND >,  G . A . KERAMTDAS 

C— — 

C— SUBROUTINE  QSOLV  FORMS  THE  SYSTEM  OF  ALGEBRAIC  EQUATIONS 


C««— ANO  CALLS  SUBROUTINE  MFART  TO  PARTITION  THE  EQUATIONS  ACCORDING 
C— — TO  THE  BOUNDARY  CONDITIONS.  IT  ALSO  CALLS  SUBROUTINES  OCOMP 


C— AND  SOLVE  FOR  THE  SOLUTION  OF  THE  SYSTEH  OF  EQUATIONS. 
C****« 

COMMON, 'BLOC 1 /  TO,  T ,  THAX, I FREQ, IBC,NUMBC, IB<2> 

COMMON/SLOC2/  Q< 222  >, YQ< 222 > 

C0MH0M/BL0C3/  H< 120), V< 1 20 >.DH< 1 20 >, 0V< 120), PR< 120> 
C0MM0N/BL0C4/  RH< 120>,RL< 1 20 >, 0X< 1 20 >, HC< 120> 

COMHOH/BLOCS/  M, WO, DT, GE,H0, XL, C2, IGEOM, IC2.AM0 
DIMENSION  GSM<  4END ) , I P< NEQ  > 

DIMENSION  ESM1 <  4, 4 ), ESM2<  4, 4  >, EF<  4  >,  NS<  4 ), BV<  2  > 

C— — 

NUHEL-N 

NPOIN-H-M 

NH9AN»< NBANO- t >/2  !  HALF  BANDWIDTH 

NORD-2  >  DEGREES  OF  FREEDOM  PER  NODE 

ND0FR-2-M0RD  !  DECREES  OF  FREEDOM  PER  ELEMENT 

4CF-NN8AN-NEQ 
4GSM-JGF-NEQ 
PI-3.14159 
C0-SQRT<  HO-GE ) 

C— — - 

DO  10  1-1,4 
DO  10  4-1,4 
ESM1( I, 4>-0. 0 
10  ESM2<  1,4  >— 0 . 0 

DO  19  1-1 ,NDOFR 
13  NS< I J-I-HORD 
C -4-4- 

DO  20  I • 1 , 4END 

20  GSH< 1 )-0 . 0 
C-4—  — 

c— —BOUNDARY  CONDITIONS 
C— — 

GO  TO  <21,22,23,24,29,26),  IBC 

21  H<  1  >— AMO 
Q<  1  )— H<  1  > 

V< HFOIN )-H<  NPOIN  >-SQRT< GE-HC<  NROIH >  >/HC<  NPOIN  > 

CKNS9  >— V<  NPOIN) 

GO  TO  30 

22  W  1  )-A«0*SlN<2.*PI-T/T0) 

0<  I  >— H<  I  ) 

V< HFOIN )— H< NFOIH >-SQRT<  GE*HC<  NPOIN ) VHC<  NPOIN  > 

0<  NEQ  >-V<  NPOIN) 

GO  TO  30 

23  V<  I  )— AMO 
9<  2  >-V<  1  ) 

N< NPOIN >-V<  NPOIN  >-MC<  NPOIN  VSQRT<r  GE*MC<  NPOIN  >  ) 
CKNEQ-1  )-N<  NPOIN) 

GO  TO  30 

24  V<1  >-AH04S!N<2.*PI-T/T0) 

0<  2  >•  V<  1  > 

M< NPOIN )-V< NPOIN )— HC< NPO IN  VSvRTf GE «MC < NPOIN  >  > 


*.*r%*V' 


0596 
0597 
0598 
0599 
0600 
0601 
0602 
0603 
6604 
0605 
0606 
0607 
0608 
0609 
061  0 
061 1 
0612 
0613 
0614 
0615 
0616 
0617 
0618 
0619 
0620 
0621 
0622 
0623 
0624 
0625 
0626 
0627 
0628 
0629 
0630 
0631 
0632 
0633 
0634 
0635 
0636 
0637 
0638 
0639 
0640 
0641 
0642 
0643 
0644 
0645 
0646 
0647 


QC  NEQ-1 >-HC  NPOIN  > 

GO  TO  30 

25  H<  1  >«AH0 
Q<  1  >-HC  1  > 

V< NPOIN >«0 . 0 
Q< NEO >-0.0 
GO  TO  30 

26  H< 1 >— AH0*SINC  2 . -PI+T/TO ) 

Q<  1  >— H<  1  > 

V< NPOIN >-0. 0 
Q< NEQ )-0. 0 
30  CONTINUE 

8V<  I  >-0< IB< 1 >> 

BV< 2 >-OC  IBC2>) 

C++*** 

C  —  —GENERATION  OF  THE  BOUNDARY  FORCE  VECTOR 
C++**- 

DO  45  LL-1 , NUMEL 
EF< 1 >-0.00 

EF< 2  >»  6 . *GE*HC LL VDXC  LL  > 

EF<  3  >— 0.  0 

EF<  4  >— 6  .  +GE+HC  LL+1  VDXC  LL  > 


40 

45 

C****+ 


NS< 1  )-NS< 1 >+NORO 
NS<  2  >— NS< 2  >*NGRD 
NSC  3 )— NSC  3  >+NORD 
NSC  4  >-NSC  4  >*NORD 
DO  40  I-t , NDOFR 

1 1 - NSC  I > 

12- JGF+II 

CSHC 1 2  >— GSHC 12  >*EFC I > 
CONTINUE 

DO  35  1-1. NDOFR 
NSC  I  )«I-NORD 


BOUNDARY  FORCE  VECTOR 


35  NSC  I )-I-NORD 
C++*** 

c—— GENERATION  OF  THE  SYSTEM  MATRICES 
C***+* 

DO  80  KK-1.NUMEL 

CALL  ELNATC  KK . ESM 1 , ESM2 > 

NSC  1  >-NSC  1  >+NORD 
NSC  2  >— NSC  2  >*NORD 
NSC  3  >— NSC  3  >+NORD 
NSC  4  >-NSC  4  >+NORD 
DO  80  I-t, NDOFR 
IX-NSC I > 

J2-JGF+II 
DO  75  J-t. NDOFR 
JJ-NSC  J> 

GSMC  J2>-CSHC  J2>*ESH2CI.J>*YQ<  JJ>  !  RIGHT  HAND  SIDE  VECTOR 
JJ- JJ-I I+NHBAH 
JS-JCSN*JJ+NEQ*II 

GSHC J5>-GSH< J5>*ESH1CI, J>  >  MATRIX  COEFFICIENT  OF  SYSTEM 
75  CONTINUE 
60  CONTINUE 


CALL  MPARTC CSHC JGSM* 1 > , GSMC JGF+t  >,NEQ,BV, IB. NBAND, NUMBC  1 
CALL  DCOMPC  GSMC  1  >,  IP.GSMC  JGSf1*1  >,  NEQ ,  NECt ,  NHBAN .  NBAND  < 

CALL  SOLVEC GSMC I >, IP, GSMC  JGF*1 >,GSMC  JGSM+1  >. NEQ , NEQ , NNBAN, NBAND > 

RETURN 

END 
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0659  C  *  **’***>#**>m>i**>*‘****— a*****:********— ******* 

0660  C***** 

0661  SUBROUTINE  XPLOT<R, T,N, IPT, IPLOT, ICPL),  N.C.CHU,  G . A . KERAMIDAS 

0662  C***** 

0663  C*****SUBROUTINE  XPLOT  PLOTS  THE  SOLUTION  ACCORDING  TO  THE  X-COORDINATE 
0664  C ***** 

0665  C0MM0N/BL0C4/'  RH<  1 20  > ,  RL<  1  20  >,  DX<  1 20  ),  HC<  1  2 0  > 

0666  C0MM0N/8L0C5/  NP, W0,DT,GE,H0,XL,CZ, IGEOM, ICZ, AMO 

0667  DIMENSION  R< N ), H< 200 >, X< 1 20 >< Y< 200  ) 

0668  C+**** 

0669  DO  10  1-1 <N 

0670  X< I)«RL< I >/1 0. 0 

0671  Y<  I  )-<H0-HC<  I  >  VI  00. 0-0.30 

0672  H< I >-R< I )/1  0 . 0-<  0 . 3-H0/1 00 . ) 

0673  10  CONTINUE 

0674  C***** 

0675  IF< IPL0T.GT.3)  GO  TO  979 

0676  IF<ICPL<GT< 0>  GO  TO  47 

0677  C+****PL0T  GRAPHICS  OUTPUT 

0678  979  WRITE  < 1 <980) 

0679  980  FORMAT< "READY  TO  PLOT  GRAPHICS,  PRESS  RETURN") 

0680  READCt,*)  IIC 

0681  IF< IGE0M.NE.3)  GO  TO  15 

0682  989  FORMAT<"^+dF%*dAf**0.110.*1 . OpO . 2q- . 3n . 25o . 05r . 01 si vC" > 

0683  WRITE< 1,981) 

0684  981  FORMAT<”l*dFt*dA(*a0.110.i»1  .  OpO .  2q~.35n.25o  .  05r  .  01  *1  vC"  ) 

0685  GO  TO  20 

0686  15  WRITE< 1.989) 

0687  20  WRITE< 1,983) 

0688  983  F0RNAT<"fc*«2h1i2j1K“) 

0689  WRITE< 1,982) 

0690  932  F0RHAT<"t*d300,0ol'St*»2a1N%*lX-AXIS%*dTM) 

0691  IF< IPT.EQ. 0)  GO  TO  <30,40,45,30.40,30,30).  IPLOT 

0692  IF<IPT<EQ<1  )  GO  TO  <47,47,47,40,45,45,40),  IPLOT 

0693  IF< IPT.EQ. 2)  GO  TO  < 47. 47 , 47 , 47, 47, 47, 45 ).  IPLOT 

0694  30  WRITE< 1,903) 

0695  GO  TO  47 

0696  40  WRITE< 1,913) 

0697  GO  TO  47 

0698  45  WR1TE< 1,923) 

0699  903  FORMAT< "%*d30, 160okS^*a2a2NL*lELEVAT10N  ^*dT" > 

0700  91  3  FORMAT<  "t*d30, 160okSL**2»2N^*l VELOCITY  > 

0701  923  FORHAT<H*d30, 160okS%**2a2N**lPRESSURE  %*dT"  > 

0702  47  WRITE< 1,984) 

0703  934  FORHAT<"l**A") 

0704  WRITE< 1,987) 

0705  987  FORMAT<"t*dF"> 

0706  C*****PLOT  COMPUTED  SOLUTION  AT  TIME  T 
0707  WRITE<  1,985)  < X<  I  >,  R<  I  ),  I- 1 ,  N  ) 

0708  WRITE<  1 , 989)  <X<H-I*1  ),H<N-I*I  >,  1-1  ,N) 

0709  IF<IPT<GT< 1)  GO  TO  48 

0710  WR!TE<  1 , 985)  <X<  I  >,Y<  I  ),  1-1  ,N> 

0711  985  FORMAT< 2F7 . 4  ) 

0712  48  WRITE< 1 ,988) 

0713  988  FORMATC "t*dE" ) 

0714  C«**** 

0715  WR1TE< 1,986) 

0716  986  FORHAT<"%*«e**dT|*dE"> 
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0717  4R I TE< 1,997) 

0719  URITE< 1,998)  T 

0719  997  FrtRMAK  "£*<1140,31  Ook  S^*m2m  1  NL*  1  ^M?.  *dTt  *m  1  m  1  N  '  > 

0720  998  F0RNAT<  *d1  40 , 31  Ook  *m2«i  t  N^*  l  T  , F5 . 3 .  M^*dTt  *m  t  nt  1  N”  > 

0721  50  RETURN 

0722  END 

0723  c*-*** 

4.KYMA  T-00003  IS  ON  CR0001  1  USING  00188  6LKS  R-0000 

0724  c****************************************» 

0725  C**-** 

0726  SUBROUTINE  TPL0T<R, TIN. N, IPT),  G . A . KERAMIDAS 

0727  c***** 

0728  C****»SUBROUTINE  TPLOT  PLOTS  THE  SOLUTION  ACCORDING  TO  THE  TIME  COORDINATE 
0729  C***** 

0730  COMMON/BLOCI /  TO. T, TMAX, IFREQ , IBC, NUMBC , IB< 2 > 

0731  C0MM0N/BL0C4/  RH< 1 20 > , RL< 1 20 >. DX< 1 20 >, HC< 1 20  ) 

0732  C0MM0N/BL0C57  N,W0,DT,GE,H0,XL,CZ, IGEOH, ICZ.AMO 

0733  DIMENSION  R<N),TIM<N) 

0734  C***** 

0735  C****'* 

0736  C.  .  PLOT  GRAPHICS  OUTPUT  . 

0737  WRITE  < 1,980) 

0738  980  FORMAT< -READY  TO  PLOT  GRAPHICS,  PRESS  RETURN") 

0739  READ<1,«)  II 

0740  WRITE< 1,981) 

0741  981  FORHAT<  ”l*dF(:*dA%*«S0. 0a3,  OptqC- ) 

0742  URITE< 1 , 982 ) 

0743  982  FORMAT< -^*d300,  0okS^*«2*i1N%*lT-AXIS^*dT- > 

0744  IF< IPT.EQ. 0)  WRITE< 1 ,903) 

0745  IF< IPT.EO. 1 )  WRITE< 1 ,913) 

0746  IF< IPT.EO. 2)  WRITE< 1 .923) 

0747  903  FORHAT< "%-d30,  160okS(+m2m2N(*lELEVATI0N  *-dT«  ) 

0748  913  FORMAT< -^*d30, 160okS^*»2n2Nt*lVEL0ClTY  L-dT* > 

0749  923  FORMAT< -t*d30. 160okS^*»2n2H^*lELEVATI0N  %-dT* > 

0750  WRITE< 1 , 984) 

0751  984  FORMATS -^*»A- ) 

0752  C.  .  PLOT  COMPUTED  SOLUTION  AT  TIME  T1  . 

0753  WRITEC1.983)  < TIH< I > . R< I > , I»1 , N > 

0754  985  FQRMAT<  2F7 . 4 ) 

0755  WRITE< 1 , 986  > 

0756  986  FORMAT< -^*«B^*dT%*dE"  ) 

0757  WR 1 TE< 1 , 998 )  T 

0758  998  FORMAT<  -%*d140,3l  Ook S^-a2m  t NL*  1  T  . F4 . 2 , “^-dT^-m 1 • 1 N“  ) 

0759  RETURN 

0760  END 

0761  c***** 


4KYMA  T-00003  IS  ON  CR00011  USING  00168  BLKS  R-0000 

•762  C*************************************** ******* 

0763  C****» 

0764  SUBROUTINE  XOUT< LUIN, IPT),  H.C.CHU 

0765  C***** 

0766  COMMON/BLOC 1 /  TO , T . TMAX . IFREQ , IBC , NUMBC , IB< 2 > 

0767  COMHON/8LOC3/’  H<  1 20  ) ,  V<  1  2 0  > ,  DH<  1 20  > ,  DV<  1  20  > ,  PP<’  1  2 0  > 

0768  COMMON^BLOC4/  RH<  1 20  > , RL< 1 20  ) , CXf 1 20  > , HC< 120) 

0769  COMMON/BLOCS/  N, WO . DT , GE , MO , XL , CZ , IGEOM, ICZ , AMO 

0770  C0MH0N/BL0C6/  HTN<  151  ),  VTN<  1 5 1  >. HT  J <  1 5 1  ) .  VTK  1 5 1  > .  T IM<  1  51  > 

0771  C***** 
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NP0IN-N*1 

IF< IPT.CT.O)  GO  TO  20 
ICPL-0 
I CPU 0-0 
WRITE< t . t  000  > 

1  000  FORMAT<  **************************** ******************< 

1  ,/,«****  DO  YOU  WANT  TO  OUTPUT  THE  RESULTS  ON  THE 

2  ,/,«****  PRINTER'*  TYPE  1  FOR  YES.  0  FOR  NO  :  “  > 
READ( LU1N. -  )  IPRIN 

IF< IPRIN.LT. t )  GO  TO  30 
URITE<  6, 1 035  > 

WRITE<6, 1003) 

1003  FORMAT<40X. -***********•*****************************« 


1  ***,/',  4 OX,  “* 
2*  -,/.40X,"* 
3*  m./.40X.m+ 
4*  " , 40X, ■* 
3*  •.✓.40X,"* 
6 *  a./.40X.m* 
7*  a , 40X, "* 
8*  *./.4QX, •* 


9*  •.✓,40X, 


COMPUTATIONAL  HYDRODYNAMICS 

SIMULATION  OF  WATER  WAVES 
IN  AN 

OPEN  CHANNEL 


C0-SQRT<N0*GE) 

CR0-C0*DT*W0 
URITE<  6, 1 060 ) 

WRITEtS, 1010)  CRO 

1010  FORMAT<40X, ************************ 
1  —  "  ,/,40X,  “**** 

2*  • 40X, ■****  THE  COURANT 

3  40X, ***** 

4*  • . /. 40X . ************************* 
3*  ”) 

WR1TE< 4,1 060 ) 

WRITE< 6, I  015)  N, XL, OX< 1 ),DT,AM0,  IBC 

1013  FORMAT< 40X, ************************ 
t  *** , / , 40X , ■**** 

. 24  .■,/,40X, ■**•*  NUMBER  OF 

3  ,/.40X, -*•**  LENGTH  OF 

4  , / , 40X, "**•*  LENGTH  OF 

3  , /, 40X, *****  TIME  STEP 

6  , ✓, 40X, *****  AMPLITUDE 

7****" , /, 40X, •***•  60UDARY  i 

92X,  ■*•**-,  A40X,  •***•  IBC 


THE  COURANT  NUMBER  IS  " , F6 . 3. 11 X . ■ < 


NUMBER  OF  ELEMENTS  IS  -,16  , 1  OX, "****■ 

LENGTH  OF  CHANNEL  IS  ■ , F6 . 3, 1 1 X, " ***** 
LENGTH  OF  ELEMENT  IS  * . F4 . 3 , 1 1 X . " ***** 
TIME  STEP  SIZE  IS  *. F6 . 3, 1 4X, ****** 
AMPLITUDE  OF  INCOMING  WAVE  IS  *  ,  F4 . 3, 2.Y.  " 
SOUDARY  CONDITION  AT  X-0.0  IS  IBC-*, 12, 
IBC-1 ,H< I )«CONSTANT*1SX. ****** 
!8C-2,H< 1 >*S1NE  WAVE  *** 

I8C-3, V< 1 >-CONSTANT  ♦*• 

IBC-4, V< 1 >-SINE  WAVE  •** 

!0C-3,M< 1  )*CONSTANT. FINITE  CHANNEL  *** 

IBC-4, M< f >-SINE  WA^C, INFINITE  CHANNEL*** 


9  , /,40X, *****  IBC-2,H< 1 >*S1NE  WAVE 

I*  *,/,40X, *****  I8C-3, V< 1  )-CONSTANT 

2*  » , /, 40X, *****  IBC-4, V< I  )«SINE  WAVE 

3*  * ,/,40X, *****  J0C-3,H< 1  )*CONSTANT. FINITE  CHANNEL 

4*  * , /, 40X, *****  IBC-4, M< I >-SINE  WA^g, INFINITE  CHANNJ 

3*  m./.40X. ***** 

4*  " , 40X,  •*•*•**•*****••***•♦♦***♦**♦***«.♦**.*♦***♦****♦***«•< 
7**"  > 

)  IF< IPRIN.LT. 1 )  GO  TO  43 
WRITE<  * . 1030)  T 
WR l  *  " .  *  1043) 

WR  .  Vi  #,1033)  <  I  ,H<  I  ),V<  I  >.PL<  I  ).PH<  I  >,HC<  I  >.PP'  I  ),  1-1  ,NPOIN> 


IF< IPT.CT.O)  GO  TO  43 
30  WRITE<  1  ,  1  020  > 

1 020  FORMAT* •****•*********•*•**#***•*••*•**♦**..♦****#♦•******♦« 

1  ,/,*•*••  DO  YOU  WANT  TO  OUTPUT  THE  RESULTS  ON  THE  TAPE 

2  ,/,•****  DRIVE'*  TYPE  1  FOR  YES.  0  FOR  NO  s  “  > 

READ< LUIN, * )  I  TAPE 


0033  43  IFC 1TAPE.LT. 1 >  CO  TO  50 

0836  UP I  TEC  8. I  04  0  >  C  HC  I  ) ,  VC  l  >,  HCC  I  ) , PRC  I  )  , RL< I > , T , 1*1 , NPOIN  > 

083?  50  IFCIPT.GT.O)  GO  TO  57 

0838  UR I TEC  1 . 1 023  ) 

0839  1  025  FORHATC  •*•*■*■■■■■••****•*»■•■*•*•■♦*■♦■**■•****■****•* 


0840 

1 

,/, •■•■■SPECIFY  UH1CH  OF  THE  FOLLOUING  VARIABLES  YOU 

*  *  *  *  “ 

0841 

4L 

-■•••WANT  PLOTTED  : 

0842 

3 

, 7, •■•■■ 

0  ■ 

NOTHING 

****** 

0843 

4 

, 7, "*••• 

1  - 

ELEVATION 

****"' 

0944 

5 

,/,-*••■ 

2  - 

VELOCITY 

****  “ 

0845 

6 

, 7, “*••• 

3  - 

PRESSURE 

****** 

0846 

7 

,7, •■*•* 

4  - 

ELEVATION  AND  VELOCITY 

**  *  *  ** 

084? 

8 

,7 , •■•■■ 

5  - 

VELOCITY  AND  PRESSURE 

****** 

09*8 

9 

,7. •*•■■ 

6  ■ 

ELEVATION  AMD  PRESSURE 

****** 

08*9 

1 

,7. "*•■« 

7  ■ 

ELEVATION, VELOCITY, AND  PRESSURE 

****  * 

0850 

2 

.7. •*•** 

IPLOT  ■  ) 

083 t  REAOCLUIN.O  IPLOT 

0832  57  IPT-0 

0853  IFC IPLOT. GT. 3)  GO  TO  58 

0834  IFC ICPL . LE . 0 .OR . ICPL0 . £Q . 1 >  CO  TO  58 

0833  UR I TEC  1.1063) 

0836  UR I TEC  1.1030) 

0837  1  030  FORHATC  *■■••••*•••••*•****■♦**•*■♦—*■*•*********••*■**•***♦***+•■" 

0858  1  .7.  •***■■  DO  YOU  UANT  A  CONTINUOUS  PLOT?  *•**“ 

0859  2  ,/,-*•••  TYPE  1  FOR  YES,  0  FOP  NO  :  *  ) 

0860  READCLUIN.*)  ICPL 

0861  ICPL 0-1 

0862  38  IFC IPLOT. LE. 0)  GO  TO  999 

0863  GO  TO  C60, 65, 70, 60. 65, 60,60).  IPLOT 

0964  60  CALL  XPLOTC H. T, NPOIN, IPT, IPLOT, ICPL  ) 

0363  IPT-IPT-M 

0366  IFC IPLOT. EO. 1 )  GO  TO  999 

0867  IFC IPLOT. EO. 6)  CO  TO  70 

0368  63  CALL  XPLOTC V . T, NPOIN . IPT , IPLOT , ICPL > 

0869  IPT-IPT-M 

0870  IFC IPLOT. EO. 4)  GO  TO  999 

0871  IFC IPLOT. EO. 2)  GO  TO  999 

0872  70  CALL  XPLOTCPR.T, NPOIN. IPT, IPLOT, ICPL) 

0873  IPT-IPTM 

0874  999  IFC ICPLO.EO. 1  )  GO  TO  100 

0875  ICPL-1 

0876  C*«*** 

0377  1033  FORHATC " 1 • > 

0978  1040  FORHATC 5X.6F 12. 6) 

0879  1043  FORHATC  777/,  •  NODAL  POINT  H-ELEVATION  V-VELOCITY 

0880  1  'COORD .  H-TOTAL  HO-INITIAL  PRESSURE  " ,  7/ > 

0881  1050  FORHATC  "  1  •. /V, 3X, 6NTIHE  -.F6.3.77) 

0882  1 053  FORHATC 1 112. 6F 1 3.6) 

0803  1060  FORHATC /7// ) 

0884  1  085  FORHATC  "l-dD*  > 

0803  1 00  RETURN 

0306  END 

0887 
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0883  ••*•••*•**•**•**•*•*•*•*+*•**•*•**•  •*•*.*■*** 

0809  C— •• 

0890  SUBROUTINE  TOUTCLUIN, NPOIN, IT, IPT ».  N.C.CNU 

0891 

0892  C0WH0N7BL0C  1  7  TO.T,  THAX,  IFREO,  I8C.NUH8C,  IBC  2  ) 

0893  C0HH0N/BL0C37  NC  120). VC  120>, DHC  12  0), DV<  120), PR <  1  20  "> 


0894  C0MM0N/BL0C4/  RH< 1 20 > , RL< 1 20 >, DX< 1 20  >, HC< 1 20 > 

0893  C0MM0H/BL0C3/  N, WO, OT, GE, HO , XL , CZ , IGEOM, IC2 , AMO 

0896  C0MM0H/BL0C6/  HTN<  151  >.VTN<  151  >,HT1<  131  ),VTK  151  >,  TIN<  151  < 

0897  C*— * 

0898  WRITE< 1 , 1040) 

0899  WRITE< 1 , 1025)  NPOIN.NPLOT 

0900  WRITE<  1 . 1 030  >  IT 

0901  WRITEi 1 , 1010) 

0902  101  0  FORNAT< 

0903  1  ,/,•——  00  YOU  WANT  PRINTER  OUTPUT  OF  THE  NUMBER  OF 

0904  2  ,/,•*—  POINTS, PLOTS, ETC.?  TYPE  1  FOR  YES,  0  FOR  NO  :  "  > 

0903  REAO<LUIN,->  JPRIN 

0906  25  IF< JPRIN, LT. 1 )  CO  TO  30 

0907  URITE<  6, 1 023  >  NPOIN,NPLOT 

0908  UR1TE<  6, 1 030  >  IT 

0909  WRITE< 6,1 033 >  XL, XL 

0910  WRITE<  6, 1 020  >  < HT1 < I >, VT1 ( I ). MTN< I », VTN< I >, TIH< I >, !• 1 , IT  ) 


091 1  30  URITE< 1 ,1015) 

0912  1013  FORMAT< -*****+*********************************+**+*+*********•*> 

0913  1  ,/,-*•**  00  YOU  WANT  TAPE  OUTPUT  OF  ELEVATION  AND  *•**• 

0914  2  ,/,■*—  VELOCITY"'  TYPE  1  FOR  YES,  0  FOR  NO  i  ") 

0913  READ< LUIN,— >  JTAPE 

0916  43  1F< JTAPE.LT. 1 )  CO  TO  50 

0917  ENOFILE  8 

0918  WRITE<8, 1020)  <  MT1  <  I  ),  VT1<  I  ).  HTN<  I  >,  VTN<  I  ) .  TIMCI  ) ,  1-1  ,  IT  > 

0919  ENOFILE  8 

0920  50  IPT-0 

0921  CALL  TPLOT<NTN, TIN, IT, IPT^ 

0922  IPT- 1 

0923  CALL  TPLOTCVTN, TIM. IT, IPT > 

0924  IPT-2 

0923  CALL  TPLOT<HTi , TIN, IT, IPT > 

0926  IPT- I 

0927  CALL  TPLOT< VT1 , TIN. IT . IPT > 

0928  C— — 

0929  1020  FORNAT<  3X , 3F 1 6 . 6 ) 

0930  1023  FORMAT<///,1 OX  .^NUMBER  OF  POINTS,  NP  •  -.I4.AV. 

0931  C  1  OX, * NUMBER  OF  PLOTS, NPLOT  •  *,I4> 

0932  1030  FORNAT< 1  OX,///, •  THE  8  OF  POINTS  FOR  THE  TINE  PLOT  IS  ••,I4.////> 

0933  1 033  FORMAU 1 3X , ■ EL E VAT I ON " . 7X , *  VEL  OC I TY ■ , 7X . " ELE VA  T I ON " . 8X , • VELOC I T Y " 

0934  1  , 1  OX. -TINE",/, 12X.-AT  X-0. 00",7X,-AT  X-0. 00", 6X. "AT  X«*. 

0933  2  F3. 1 ,6X, -AT  X--.F3. 1//) 

0936  1040  FORHAT<  ) 

0937  C— — 

0938  RETURN 

0939  END 
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£*+**• 

(»«*•*« 

c*«— 


c— — 
C— — 
c*-*— 
c— — 

c*—*« 
c— — — 
c»*«— 


SUBROUTINE  PARTH<  CSN.  CP.  NP,  BV.  IB.  N8W.  NBC  > .  G  .  A  . KERAMIDAS 

SUBROUTINE  FOR  PARTITIONING  OF  THE  GLOBAL  STIFFNESS 
MATRIX  t  FORCE  VECTOR  BY  THE  METHOD  OF  DELETION  OF 
ROWS  6  COLUMNS  ACCORDING  TO  THE  SPECIFIED 
BOUNDARY  CONDITIONS. 

THE  GLOBAL  STIFFNESS  MATRIX  G$H<  NP , N6AND  >  CONTAINS 
THE  NBANO  NON-ZERO  DIAGONALS  OF  THE  < NP-NP >  MATRIX. 
<<WARNING>>  PAPTM  IS  USED  ONLY  FOR  ( EMA  >  OPRAYS  • 
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0019 

DIMENSION  9V<  N9C  >, I9< N9C  >, ND( 1?  > 

0010 

EHA  GSM<  NP , N9W  > ,  GF  <  HP , 1  ) 

OOt? 

c*— • 

0019 

N3«<N9W-1  i/2 

0019 

NP4-NP-2-N9C 

0020 

c**«— 

0021 

DO  10  I  —  1  , N9M 

0022 

10 

ND<  1  >-I 

0023 

c*— • 

0024 

ID-0 

0023 

DO  200  L-1.N9C 

002* 

ID-ID-1 

002? 

IP-I9<L> 

0029 

9C-9V< L > 

0029 

c— — 

0030 

MODIFICATION  OF  TNC  GLOBAL  3TIFFNESS  MATRIX 

0031 

c«— — 

0  THE  FORCE  VECTOR 

0032 

c— 

0033 

NS-N9U-N3 

0034 

MP-IP-NS 

0033 

IF<HF.CT.NF>  GO  TO  100 

0030 

00  ISO  1-1 / NS 

0037 

IF< I . EQ . IP  )  GO  TO  130 

0030 

ji-ip-i-ns 

0039 

CF< 1,1 >— CF< I , 1 >-CSH< I , J1 >*GF< IF, 1 >/GSM< IP,NS> 

0040 

DO  120  J-I,NS 

0041 

IF<  J . EQ . IF  >  GO  TO  120 

0042 

JJ-J-IP-NS 

0043 

I J-J-I-NS 

0044 

G3H< 1 ,  I  J >-GSM< I ,  I  J >-GSM< f  ,  J 1  >*GSM< IP, JJ>/GSN< IP, NS) 

0043 

120 

CONTINUE 

0040 

130 

CONTINUE 

0047 

CO  TO  193 

0049 

100 

DO  190  1-1, NS 

0049 

IN— I—NP4 

0030 

IF< IN.EO. IF  >  GO  TO  190 

0031 

JL-IP— IN—NS 

0032 

GF< IN, 1 )-GF< IN , 1 )-GSH< IN , JL)-GF<1P.  1 VGSM< IP. MS'* 

0033 

DO  170  J— 1 , NS 

0034 

JN- J—NP4 

0033 

IF< JN.E0.IP>  GO  TO  170 

0030 

IL-JN-IP-NS 

003? 

JC-JN-IN-NS 

9039 

GSH< IN,  JC>-GSH< IN, JC  )-GSM< IN, JL >-G3H< IP, IL >/GSM< IP, NS 

0039 

170 

CONTINUE 

0900 

190 

CONTINUE 

0001 

193 

GF< IP, 1 >-GSN< IF,N3>-9C 

0002 

NP0-IP-N8 

0903 

IF<NP0.GT.NP)  GO  TO  193 

0904 

DO  190  10-1 ,NS 

0003 

IF< 1 0 . CO . IF  >  GO  TO  190 

9900 

11-IP-IO-NS 

0907 

I2-10-IP-NS 

0909 

GSN< IF, 12 >-0.0 

0009 

com  lo.n  >-o.o 

0970 

190 

CONTINUE 

0071 

CO  TO  200 

9972 

193 

DO  199  JO— 1 , NS 

0073 

JN0-J9-NP4 

9974 

IF< JNO.EO. IP)  CO  TO  190 

9973 

J1-IP-JN0-N* 

0070 

J2- JN 0-IP-NS 

9977 

C3H< 17. J2 >-0. 0 

0079 

GSN<  JNO, J  t j«0 . 0 

0879 

196 

CONTINUE 

0080 

280 

CONTINUE 

089> 

RETURN 

0892 

END 
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8882 

8803 

9U9R0UTINE  DCON7<  9. I NT, A. N. MAXN. Ht , M3 ) ,  A  . J  . BAKER 

0994 

c ***** 

0909 

CONFUTE  THE  LU  DEC0N70S I T I ON  OF  A  NONSYNNETRIC 

0009 

BANOED  MATRIX  STORED  IN  DIAGONAL  FORM  BY 

9997 

USING  7ARTIAL  71 VOTING. 

9999 

9999 

DIMENSION  9< NAXN . M 1  >.  A<  MAXN.  M3  > .  INK MAXN  > 

9919 

901 1 

L»H1 

0012 

DO  129  I«1 , Ml 

0013 

K2»M1+2-I 

9914 

DO  100  J-K2.H3 

9913 

100 

A<  I ,  J-L  >-A<  I ,  J  ) 

9019 

L-L-l 

9017 

K2»H3-L 

0019 

DO  110  J-K2.N3 

0919 

110 

A<  I .  J)»0. 0 

0920 

120 

CONTINUE 

0021 

L«N1 

0022 

DO  220  K«1 ,M 

9923 

X»A<  K , 1 > 

9924 

1«K 

0023 

K2«K41 

9029 

IF<L.LT.N>  L-L-M 

9027 

I7<  L . LT , K2 >  GO  TO  150 

9029 

DO  149  J-K2.L 

0029 

IF<  A98< A<  J , 1 )  )-A9S< X ) )  140.140.130 

9939 

130 

X«A< 4, 1 > 

0931 

1«J 

0032 

140 

CONTINUE 

9933 

130 

INT<  K  >«I 

9934 

IF<  X  >  199.230.100 

0933 

190 

I7<!~K>  170. 190, 170 

9939 

170 

DO  199  J-1.N3 

9937 

X«A<K, J) 

9939 

A<K.4)*A<  I.  4> 

9939 

190 

A<  I ,  J  >-X 

0040 

199 

1FCL.LT ,K2>  CO  TO  220 

9841 

DO  219  J*K2.L 

9942 

X»t .  0 

9043 

IF<A<K,t>.NE.0.0>  X»A<  J.  1  VA<K.  1  » 

9944 

9<  K ,  J-K  >*X 

9943 

DO  209  44»2.ni 

9049 

200 

A<  4 ,  44-t  >«A<  J.  44  >-A<  K.  44  >** 

9947 

210 

A<  J,H3>«0.0 

9949 

220 

CONTINUE 

9049 

RETURN 

9939 

9931 

A  EERO  7IV0T  HAS  BEEN  FOUNO 

9932 

7R1NT  ERROR  MESSAGE  AND  ST07 

9933 

(**«•« 

9934 

230 

CONTINUE 
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0035  WRITE<  6, 600  >  K,  I ,  K2 ,  L .  J,  A<  K,  1  <,*•.  j,  1  >.:< 

0036  600  FORMAT</V,  ">>>>>  ZERO  PIVOT  IN  DCOMP  <*<»<",//.  . 

0037  t  I  OX,-  K  ,  I  .  K2  ,  L  J  ,  m<K,I>  ,  AO.I) 

0033  2  X-.//,  10X,3I5,3F10.4) 

0033  STOP 

0060  END 
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0002  C-— * 

0003  C— —*——•*— •—•———*****— —****—— —***—****• 
0004  C*+*** 

0003  SUBROUTINE  SOLVE<S,  IHT.F,  A,  N, MAXN, Ml  ,H3>,  A  .  0  .  BAKER 

0006  C— •— 

0007  C— —  COMPUTE  THE  SOLUTION  OP  THE  LINEAR  SYSTEM 
0008  C— —  USING  THE  LU  DECOMPOSITION  OP  THE  STIFPNESS 
0009  C— —  MATRIX  PROVIDED  BY  DCOMP  AND  THE  RIGHT 
0010  C— —  HAND  SIDE  VECTOR. 

0011  C*—— 

0012  DIMENSION  S< MAXN, Ml > , F< MAXN >, A< MANN. M3 > , INK MAXN > 

0013  C-— — 

00*4  L-M1 

0013  00  130  K-I.N 

0016  I-INT<  K  ) 

0017  IF<  I-K»  1  00,  1  1  0,  100 

0018  100  X-P<K> 

0019  P<K)-P<I) 

0020  P< I )— X 

0021  110  K2-K-M 

0022  IF<L.LT.M>  L-LM 

0023  IP<L.LT.K2>  GO  TO  130 

0024  DO  120  I-K2A 

0023  X-S<K,I-K> 

0026  120  F< I »-P< I >-X*F<  X  > 

0027  130  CONTINUE 

0028  L-t 

0029  DO  180  I 1—1 , N 

0030  I -N4  1  - 1 1 

0031  X-P<I> 

0032  M— 1-1 

0033  IF<  L-1 >  140,160,140 

0034  140  DO  ISO  K-2,L 

0033  130  X-X-A<  I . X  >*F<  K*M  > 

0036  160  F< I )-X/A< 1,1) 

0037  IP<L-M3)  170,180,180 

0038  170  L-L*1 

0039  180  CONTINUE 

0040  RETURN 

0041  END 
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